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ABSTRACT 

Current  numerical  barotropic  and  baroclinic  forecasting  models  have 
shown  an  inability  to  correctly  represent  the  behavior  of  the  planetary- 
scale  waves.   A  mathematical  model,  based  on  the  thermodynamic  energy 
equation,  has  been  proposed  by  Wiin-Nielsen  as  a  means  of  predicting  the 
behavior  of  the  planetary-scale  waves.   The  model  is  programmed  for  the 
Control  Data  Corporation  1604  digital  computer  to  utilize  the  operational 
output  of  machine-processed  data  and  analyses  produced  by  the  U.  S.  Navy 
Fleet  Numerical  Weather  Facility  (FNWF) ,  Monterey,  California.   The 
present  form  of  the'  program  yields  a  1-hr  planetary-scale  forecast  at 
1000  mb,  850  mb,  700  mb,  500  mb  and  300  mb  in  ten  minutes. 
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1.    Introduction. 

Current  numerical  barotropic  and  baroclinic  forecasts  of  the  height 
of  the  500-mb  surface  show  an  inability  to  correctly  represent  the 
behavior  of  planetary-scale  waves.   Early  models  caused  these  long  waves 
to  retrogress  at  speeds  approaching  the  Rossby  speed.   This  phenomenon 
is  not  observed  in  the  atmosphere  and  much  work  has  been  done  in  an 
effort  to  correct  this  difficulty.  Wolff  [7]  has  shown  that  forecasts 
by  the  non-divergent  barotropic  model  are  more  nearly  correct  if  the 
planetary-scale  waves  are  held  stationary.   Another  approach  to  correct 
this  defect  has  been  to  introduce  a  Helmholtz  term  into  the  vorticity 
equation.   This  is  equivalent  to  adding  a  small  amount  of  divergence. 
Wiin-Nielsen  [5]  pointed  out  that  an  analysis  of  baroclinic  instability 
and  phase  speeds  of  waves  in  a  baroclinic  atmosphere  leads  to  the  con- 
clusion that  the  phase  speed  for  large  wave  lengths  approaches  the 
Rossby  speed  in  a  non-divergent  atmosphere. 

Burger  [l]  has  shown,  from  scale  considerations,  that  the  vorticity 
equation  loses  its  prognostic  value  and  assumes  a  quasi-stationary 
character  for  waves  whose  lengths  are  of  the  same  order  of  magnitude  as 
the  radius  of  the  earth.   Adopting  this  as  a  working  hypothesis,  and 
defining  waves  of  this  magnitude  as  planetary-scale  waves,  Wiin-Nielsen 
[6]  developed  a  prediction  model  for  the  planetary-scale  waves. 

This  investigation  is  an  attempt  to  provide  a  means  of  correctly 
predicting  the  planetary-scale  motion  based  on  the  Wiin-Nielsen  model. 
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2.   The  Model. 

The  vorticity  equation  for  large-scale  atmospheric  flow  in  isobaric 
coordinates  is 

Burger  [l]  showed  that  this  equation  may  be  greatly  simplified  by  order- 
of -magnitude  considerations  when  applied  to  the  planetary  scale. 
Equation  (1),  after  Burger,  reduces  to  the  stationary  form 


/2sr  ^"yv-v-o . 


(2) 

Wiin-Nielsen  [6]  demonstrated  that  equation  (2)  leads  to  instability  on 
all  scales,  with  large  amplification  rates  for  the  smaller  scale  motions 
which  may  be  present.   The  shorter  waves  were  stabilized  by  adding  the 
advection  of  relative  vorticity  into  equation  (2) ,  leading  to 


V'V/i  ^S^V'V-  0 


or  (3) 

The  coriolis  parameter  in  equation  (2)  is  modified  to  +    in  equation  (3) 
to  avoid  mean -vorticity  generation  at  any  level.   The  inclusion  of  rela- 
tive vorticity  leads  to  the  delineation  of  a  region  of  instability.   In 
the  meteorological  range  of  vertical  wind  shear  and  wave  length,  the 
shortest  e-folding  time  is  of  the  order  of  5  days.   The  e-folding  time 
is  the  time  required  for  the  amplitude  to  increase  by  a  factor  of  Q, 
Since  the  amplification  rate  is  given  by  ^  »—    ,   the  e-folding  time 
is  simply 


The  thermodynamic  energy  equation  is  derived  in  a  straightforward 
manner  from  the  first  law  of  thermodynamics ,  applied  to  adiabatic  flow^ 
giving 


ja-fVv-vQfV^[p|Cij-^Ii 


=  0.  (4) 


This  equation  may  be  transformed  to  equation  (5)  by  approximating  sLs 
by  finite  differences. 


where 

cr   = 


To  apply  equations  (3)  and  (5),  the  atmosphere  is  divided  into 
9  layers  by  isobaric  levels  numbered  0  to  9,  where  the  levels  0  and  9 
correspond  to  200  mb  and  1000  mb,  respectively  (see  figure  1) , 

mb     Input     Interpolated     Computed     Forecast  k 

200    —  [) 0 

300    —  D  — -  Yi 1 

400   —  [)  =--« <r,  J:^ ^  .._  2 

500    —  D  Y[ 3 

600  —  D ^"-^   —  k  — -  4 

700    —  D  71  ^~~ •  5 

775    1^ C:o> 1^ 6 

850    — D 71    — 7 

925    —  Q  — xr,i^   k  8 

1000    —  [) — —  2.  9 

Figure  1.   Level  arrangement  for  the  model. 

-3- 


The  vertical  velocity  is  computed  at  even  levels  by  applying  the 
vorticity  equation,  (3),  at  the  odd  levels  as  follows; 

Here  k  is  even  and  P   is  the  pressure  difference  between  even  levels. 
The  vertical  velocity  is  assumed  to  be  zero  at  200  mb. 

The  prognostic  equations  are  obtained  by  applying  the  thermodynamic 
energy  equation,  (5),  at  the  even  levels,  namely, 


l>— ^.,-Vk^H 


5^5" 


Equation  (6)  may  also  be  expressed  as 


■=r=   tlU-,] 


(8) 


•-    ■     r^-T-  KU"  0)  (9) 


or 


Equation  (9)  is  used  to  evaluate  the  stability  parameter  in  a  standard 
atmosphere  at  the  even  levels.   These  values  are  included  in  the  program 
as  constants. 

The  lower  boundary  condition  is  established  by  noting  that 

The  second  term  of  the  right  hand  side  of  equation  (11)  is  zero  in 


geostrophic  flow.   Since  W    at  the  surface  is  given  by 

equation  (11)  may  be  rearranged  to  obtain  equation  (12a) j 

^t         ^r     V^    ^   "^  (12a) 

or,  approximately, 

>J^  i-  Itt    ^   -L^el^*  ^^-VHt  ,  (12b) 

where  the  surface  wind,  vy    ,  has  been  approximated  by  V/a        .      The 
difficulties  embodied  in  equations  (12a)  and  (12b)  are  rather  serious 
since  i-0   at  the  lower  boundary  is  usually  obtained  by  neglecting  ^-f 
in  equation  (12a).   Thus,  it  is  likely  that   £3     is  a  small  differ  =■ 
ence  of  two  relatively  large  terms.   In  spite  of  these  misgivings ,  the 
initial  program  was  designed  in  accordance  with  Wiin-Nielsen°s  model  as 
approximated  by  equation  (12b) . 

The  computations  carried  out  during  this  investigation  were  made 
with  a  Control  Data  Corporation  1604  (CDC  1604)  digital  computer,  with  a 
core  storage  capacity  of  32,768  words  of  48  bits  each.   The  input  data 
were  obtained  from  the  bi-daily  analysis  of  the  height  of  the  300-mb, 
500-mb,  700-mb,  850-mb  and  1000-mb  surfaces  by  the  U.  S.  Navy  Fleet 
Numerical  Weather  Facility  (FNWF) ,  Monterey,  California. 

In  order  to  utilize  available  data  and  computer  subroutines  from 
FNWF,  the  equations  were  converted  to  finite-difference  form  and  scaled 
for  fixed-point  fractional  arithmetic. 

A  1977-point,  square-net,  octagonal  grid,  circumscribed  by  latitude 

ION  on  a  polar  stereographic  projection  of  the  Northern  Hemisphere,  was 

employed  to  represent  the  data  fields.   The  net  spacing  is  381  km  at  60N 

where  the  projection  is  true. 
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The  scaled,  finite-difference  equations  corresponding  to  equations 
(8),  (7),  (10)  and  (12)  are  equations  (13),  (14),  (15)  and  (16), 
respectively. 


^ 


5:  --  t'[(iHo>'(?,  f^MK      a., 


The  scaling  conventions  used  in  these  equations  are 


lO  = 

u3.  Z 

p- 

D- 

P  ' 

?  -2- 

m  ■- 

A         ^' 

n^ 

^■z^ 

7  = 

i-£' 

cr  - 

^•2" 

(tnb  -  sec   ) 


(sec  ) 


(cm) 


cm  sec   mb   ) 


(16) 


(cm) 

(mb)  (17) 

(non-dimensional) 


(sec  ) 


K  ■-  k  'C  (cm) 

Equation  (6)  represents  the  standard  stability  plus  an  anomaly  due 

A 

to  the  departure  from  the  standard  atmosphere.   The  quantity   O""   in 
equation  (15)  refers  to  the  stability  parameter  at  level  k  in  a  standard 
atmosphere.   The  remainder  of  the  right  hand  side  of  equation  (15)  is  the 
finite-difference  form  of  equation  (10) j  which  was  used  to  approximate 
the  stability  anomaly  at  each  grid  point.   The  point  values  of  stability 
are  averaged  over  each  pressure  surface  to  maintain  consistency  with 
respect  to  the  energetics  of  the  model. 

The  vertical  derivatives  in  the  calculation  of  the  stability  anomaly 
are  approximated  in  a  conventional  manner  [3s,  pp.  47-49].   The  height 
fields  at  the  even  levels  (figure  1)  are  approximated  by  a  simple  linear 
relationship.   It  is  assumed  that  the  ratio  of  the  height  of  a  level  to 
its  adjacent  level  in  a  standard  atmosphere  is  the  same  as  the  ratio  of 
the  heights  of  these  levels  in  the  real  atmosphere.   A  test  of  the 
validity  of  this  assumption  indicates  that  the  standard  deviation  of  the 
actual  ratios  with  respect  to  the  standard  ratios  is  two  orders  of  magni- 
tude less  than  the  ratios  themselves.   The  approximation  leads  to  a  maxi- 
mum error  of  97o  of  the  height  in  the  54  cases  examined^  with  an  average 
error  of  less  than  1%.   The  approximation  is  less  accurate  in  the  lower 
portions  of  the  atmosphere  than  it  is  in  the  upper  regions  below  the 
tropopause. 

Additional  boundary  conditions  are  imposed  on  the  finite-difference 
approximation  to  the  model.   These  include- 

a)   setting  the  forecast  height  and  thickness  changes  at  the  grid 

boundary  points  to  zero; 
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b)   setting  the  relative  vorticity  and  vertical  velocity  at  the 

grid  boundary  points  to  zero. 
The  order  of  computation  during  a  time  step  is  as  follows" 

a)  the  previously  prepared  planetary-scale  height  analyses  at  the 
five  odd  levels  (see  figure  1)  are  entered  into  the  computer 
and  the  height  values  at  the  five  even  levels  are  approximated; 

b)  the  stability  parameter  is  computed  at  each  grid  point  on 
levels  2,  4,  6  and  8  and  is  pressure  averaged  for  each  level; 

c)  the  vertical  velocity  is  computed  at  each  even  level  from  the 
advection  of  absolute  vorticity  at  the  next  lower  odd  level 
after  establishing  the  vertical  velocity  at  level  0  as  zero; 

d)  the  thickness  fields  between  adjacent  odd  levels  are  obtained 
from  the  input  height  fields  and  the  advection  of  thickness  is 
calculated; 

e)  the  advection  of  thickness  and  the  vertical  velocity  term  at 
each  even  level,  except  level  Oj  is  combined  to  form  the  local 
thickness  change; 

f)  the  height  change  at  level  9  is  determined; 

g)  the  calculated  height  and  thickness  changes  are  combined  with 
the  initial  height  at  level  9  and  the  thicknesses  of  the  upper 
layers  to  form  the  forecast  height  fields  at  levels  1,  3$  5g 

7  and  9. 
Forward-difference  approximations  are  used  in  the  initial  time  step. 
All  subsequent  time  steps  employ  central  differences. 
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3.    Investigative  Procedures  and  Results o 

An  initial  attempt  at  evaluating  the  model  was  made  by  providing 
the  space-mean  height  fields  at  300  mbj  500  mb^  700  mb;,  850  mb  and 
1000  mb  as  input.   These  fields  are  obtained  by  passing  the  analyzed 
height  fields  through  four  scans  of  the  smoothing  equation, 

0^  J^^  .  r-^-^    A^"^ 


(18) 


In  this  instance,  k  =  1/8,   This  approach  proved  unsuccessful.   An 
analysis  by  Holl  and  Clark  [2]  of  the  wave -component  residuals  shows 
that,  after  smoothingj  507o  of  the  energy  of  wave  number  16  remains  in 
the  space-mean  pattern  at  latitude  45N,.   This  is  sufficient  justifica- 
tion to  remove  the  FNWF  space-mean,  in  its  present  form;,  as  a  suitable 
input  for  the  model  under  consideration. 

The  1000-mb  height-change  forecasts  produced  by  the  approximations 
of  equations  (12a)  and  (12b)  were  unsatisfactory.   The  forecast  changes 
had  no  pattern  and  were  an  order  of  magnitude  larger  than  observed 
values.   This  result  casts  further  doubt  on  the  utility  of  equations 
(12a)  and  (12b). 

A  second  approach  was  to  provide  the  model  with  input  fields  which 
were  prepared  by  a  hemispheric  Fourier  analysis  of  the  height  fields. 
The  zonal  mean  and  the  components  due  to  waves  1^  2  and  3  were  extracted 
and  added  together.   These  summed  fields  provided  a  planetary-scale 
input  to  the  model.   Figure  2  is  a  sample  of  one  of  the  fields  used  in 
the  limited  evaluation  of  the  model. 

The  first  attempts  with  the  Fourier  input  were  not  successful. 
Short-wave  components  were  generated  and  first  appeared  in  the  mapped 
height  fields  after  a  2-hr  forecast.   These  components  continued  to  grow 


in  amplitude  and  progressed  around  latitude  circles,  as  shown  in  figure  3 
for  an  S-hr  forecast „   These  short-wave  phenomena  are  tentatively 
attributed  to  computational  instability.   The  time  interval  used  in 
these  unstable  forecasts  was  one  hour„ 

The  1000-mb  height-change  forecast  was  not  unstable;  however ^  the 
earlier  difficulties  persisted.   On  this  basis^  the  utility  of  equations 
(12a)  and  (12b)  was  rejected.   In  the  interest  of  evaluating  the  main 
features  of  the  models  a  simple  approximation  to  the  lower  boundary 
conditions  was  adopted.   The  12-hr  height  changes  at  the  1000-mb  surface 
were  determined  from  the  result  of  the  Fourier  analysis  of  the  initial 
height  fields.   It  was  then  assumed  that  the  appropriate  fraction  of  this 
change  would  be  applied  to  the  1000-mb  surface  at  each  time  step   Such 
an  approximation  leads  to  minor  inconsistencies  during  the  iterative 
procedure,  but  leads  to  a  perfect  12-hr  forecast  at  1000  mb. 

The  problem  of  possible  computational  instability  has  been 
approached  on  the  basis  of  the  general  considerations  outlined  by 
Thompson  [4],   The  1-hr  time  interval  was  halvedj  while  retaining  the 
approximation  for  the  lower  boundary. 

At  the  time  of  writing,  26  time  steps  have  been  calculated  with  no 
sign  of  short-wave  phenomena,  which  appears  to  support  the  earlier  sug= 
gestion  of  computational  instability  associated  with  1-hr  time  steps. 
Figure  4  is  a  12-hr  forecast  for  the  500-rab  level  from  12Z  19  Novem- 
ber 1962,   The  errors  associated  with  this  forecast  are  shown  in 
figure  5,   The  forecast  12-hr  height -change  field  is  figure  6  and  the 
actual  12-hr  height-change  field  is  figure  7,   The  analyzed  500-mb  field 
for  OOZ  20  November  1962  is  figure  8,   The  pillow  and  root -mean-square 
error  from  figure  6  are  very  similar  to  those  obtained  from  figure  7» 
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(Pillow,  P5  and  root -mean-square  error,  RMSEj  are  defined  in  Appendix  I.) 
In  figure  6,  P  =  -2  feet  and  RMSE  s  71  feet<,  while  figure  7  gives 
P  s  6  feet  and  RMSE  s  74  feeto   An  examination  of  figure  5  shows  a  region 
of  large  errors  in  the  vicinity  of  the  North  Pole  (at  the  center  of  the 
figure) ,   This  is  to  be  expected  since  the  features  of  the  imput  data  are 
not  of  planetary  scale  at  these  latitudes.   The  major  errorSj,  with  a 
maximum  of  -480  feet,  are  associated  with  the  slightly  deepening  low 
center  in  the  northwest  region  of  figure  2.   The  model  forecasts  a 
lowering  of  central  height  to  16,560  feet.   In  actuality,  the  height 
lowered  to  16,770  feet  at  the  position  one  grid  point  to  the  northeast,. 
Most  of  the  other  errors  in  the  lower  latitude  regions  are  less  than 
100  feet.   These  errors  appear  small  enough  to  warrant  further  evaluation 
of  the  model.   Appendix  I  is  a  tabulation  of  the  pillow  and  RMSE  of  error 
fields  generated  when  osing  persistence  as  a  forecast  on  selected  days  of 
November  1962.   These  indicate  the  forecast  improvement  required  of  the 
Wiin-Nielsen  model. 
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4,   Suggestions  for  Future  Study. 

The  limited  results  described  in  section  3  are  not  yet  adequate  to 
evaluate  the  usefulness  of  the  long-wave  models   There  remain  several 
tasks  to  be  accomplished  before  this  may  be  established. 

a)  Revise  and  streamline  the  present  computer  program  to  reduce 
the  amount  of  time  required  to  make  the  calculations  for  one 
time  step.  At  present,  this  time  is  5  minutes;  however ^  this 
may  be  reduced,  by  more  sophisticated  tape  handlings  to  about 
3  minutes.   Further  reduction  might  be  made  by  packing  data  in 
the  computer  memory  to  minimize  or  eliminate  the  necessity  of 
data  transfer  to  and  from  magnetic  tape. 

b)  Conduct  an  extensive  evaluation  of  the  model  on  many  sets  of 
data  to  prove  or  disprove  the  adequacy  of  equations  (3)  and 
(5). 

c)  If  the  evaluations  of  b)  prove  the  model  to  be  successful s  a 
suitable  lower  boundary  condition  must  be  devised  and  tested 
in  order  to  obtain  operational  planetary-scale  forecasts.  A 
further  desirable  feature  for  an  operational  forecasting  model 
would  be  the  elimination  of  the  required  Fourier  analysis  of 
the  input  data.  The  most  likely  possibility  here  may  be  the 
use  of  one  of  the  new  smoothing  techniques  presently  undergo 
test  at  FNWF, 
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Figure  2  -  500-mb  height  analysis  of  the  Northern 
Hemisphere  after  Fourier  analysis, 
12Z  19  November  1962. 
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Figure  3  - 


500-mb  planetary-scale  height  forecast 
for  the  Northern  Hemisphere  showing 
features  e^valuated  aa  computational 
Instability,  8-hr  forecast  from  12Z 
19  November  1962.  1-hr  Increments. 


Figure  4  -  12-.hr  500-mb  planetary-scale  height  forecast 
for  the  Northern  Hemisphere  from  12Z   19 
November  1962.  30-min  increments. 
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Figure   5 


12-hr   500-mb   planetary-scale   height -forecast 
errors    for   the   Northern  Hemisphere    from 
122   19  November    1962.    30-mln   Increments. 


Figure  6  - 


12-hr  500-mb  planetary-scale  height -change 
forecast  for  the  Northern  Hemisphere  from 
12Z  19  November  1962.  30-min  increments. 


Figure  7     - 


12-hr   500-mb   planetary-scale    height-change 
analysis    of   the   Northern  Hemisphere    from 
12Z   19  November    1962. 
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Flgiirer  8  -  500-mb  height  analysis  of  the  Northern 
Hemisphere  after  Fourier  analysis » 
OOZ  20  November  1962. 


Figure   9     -      Polar   atereographlc   pro.lectlon 
of  the   Northern  Hemisphere   to 
the   scale   of  figa.    2  through  80 
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SCALE:  1:120,000,000 
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APPENDIX  I 


VERIFICATION  OF  PERSISTENCE  FORECASTS 
FOR  SELECTED  DAYS  OF  NOVEMBER  1962 


Pillow  is  defined  as 


r  =  -^ .  Ci'i) 

The  root-mean-square  error  (RMSE)  is  determined  as  follows: 


RMSE  = 


>-6)-P3, 

(1-2) 


^L(A-8)-P], 


In  equations  (I-l)  and  (1-2)5  A   and   Q   represent  the  point 
values  of  the  fields  being  compared.   In  this  examples  when  preparing 
the  data  for  these  tables,  the  A-field  was  the  analysis  at  the  initial 
time  and  the  B-field  was  at  a  subsequent  time. 


Index  of  Tables 
Table  I    -   Pillow  and  RMSE  for  12-hr  forecasts 
Table  II   -  Pillow  and  RMSE  for  24-hr  forecasts 
Table  III  -  Pillow  and  RMSE  for  36^  48,  60  and  72-hr  forecasts 

Table  format   -   Pillow  /  RMSE  (values  in  feet) 
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Table  I.   Pillow  and  RMSE  of  persistence  verifications  at  5  levels 
for  a  12-hr  forecast  during  November  1962 


Initial  Time 

Level 

300mb 

500mb 

700mb 

850mb 

lOOOmb 

12Z  19 

22/99 

6/74 

0/57 

-2/50 

0/61 

OOZ  20 

-23/97 

-9/68 

3/54 

2/57 

-12/69 

12Z  20 

-11/88 

-3/71 

-3/63 

0/66 

10/74 

OOZ  21 

1/96 

5/74 

4/51 

4/45 

7/50 

12Z  22 

4/87 

5/61 

5/48 

7/42 

6/50 

12Z  23 

0/90 

5/76 

-6/73 

-3/71 

5/68 

OOZ  24 

-5/95 

0/65 

6/63 

5/69 

0/74 

12Z  24 

-3/97 

-4/77 

-7/71 

-2/70 

4/66 

Table  II.   Pillow  and  RMSE  of  persistence  verifications  at  5  levels 
for  a  24-hr  forecast  during  November  1962. 


Initial  Time 

Level 

300mb 

500mb 

700mb 

850mb 

lOOOmb 

12Z  19 

-1/143 

-3/96 

2/77 

0/89 

"11/106 

OOZ  20 

-35/122 

-13/95 

0/92 

2/106 

-1/116 

12Z  21 

-16/135 

1/100 

7/90 

12/94 

18/89 

12Z  22 

0/122 

10/91 

14/82 

13/83 

14/91 
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Table  III.   Pillow  and  RMSE  of  persistence  verifications  at  5  levels 
for  36,  48,  60  and  72-hr  forecasts  during  November  1962 


Initial 

Forecast 

Time 

Interval 

Level 

300mb 

500mb 

700mb 

850mb 

lOOOmb 

12Z  19 

36 

-13/159 

-6/115 

0/116 

0/130 

0/155 

12Z  22 

36 

0/139 

5/140 

7/134 

10/133 

20/139 

12Z  19 

48 

-11/188 

-1/147 

4/137 

4/150 

6/171 

12Z  22 

48 

-4/144 

6/151 

14/163 

16/174 

20/185 

I2Z  22 

60 

-8/178 

1/188 

7/203 

13/208 

24/213 

12Z  22 

72 

-17/231 

4/223 

20/223 

21/227 

23/233 
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